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In this paper we present two new procedures for automatic detection and picking of P-wave 
arrivals. The first involves the application of kurtosis and skewness on the vector magnitude of three 
component seismograms. Customarily, P-wave arrival detection techniques use vertical component 
seismogram which is appropriate only for teleseismic events. The inherent weakness of those methods 
stems from the fact that the energy from P-wave is distributed among horizontal and vertical 
recording channels. Our procedure, however, uses the vector magnitude which accommodates all 
components. The results show that this procedure would be useful for detecting/picking of P-arrivals 
from local and regional earthquakes and man-made explosions. The second procedure introduces a 
new method called "Ratios in Higher Order Statistics (RHOS)." Unlike commonly used techniques 
that involve derivatives, this technique employs ratios of adjacent kurtosis and skewness values to 
improve the accuracy of the detection of the P onset. RHOS can be applied independently on 
vertical component seismogram as well as the vector magnitude for improved detection of P-wave 
arrivals. 

PACS numbers: 91.30.-f, 93.85.Rt 



I. INTRODUCTION 



Detection of P-wave arrivals accurately is an important 
step to make earthquake forecast, to analyze the interior 
structure of the earth and to study earthquake sources. 
Picking seismic phase arrivals correctly is also helpful 
to discriminate between natural earthquakes and man- 
made explosions. Several P-wave detection techniques 
utilize only one component of the usually recorded three- 
component seismogram [l|. These techniques work very 
well for the case of teleseismic events (with epicentral dis- 
tance A > 3300 km) since seismic P-waves approach the 
seismic sensors in a nearly vertical incidence (where A is 
the epicentral distance between a seismic station and the 
seismic epicenter) (Figure 2). 

Magotra 2] attempted to use horizontal component 
seismograms in order to determine the direction of ar- 
rival of seismic waves towards a seismic station. l3| used 
kurtosis and skewness of the vertical component seismo- 
gram to estimate the P-onset time assuming that the 
maximum of the derivatives of the kurtosis and skew- 
ness curve, just before the curve attains its maximum, 
is considered to be the time of arrival for the P-wave. 
Unlike in teleseismic events, the P-wave signal originat- 
ing from regional (100 km < A < 1400 km) and local 
distances A < 200 km) arrive at the seismic station with 
significant strengths in both horizontal and vertical di- 
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rections 0-[Jl. Since the energy is distributed between 
the horizontal and the vertical recording channels, both 
horizontal and vertical component seismograms need to 
be taken into account for improved detection. In this 
study the seismic vector magnitude of three component 
seismograms from a single station is used to determine 
the P-wave arrival without employing the derivatives of 
the kurtosis and skewness. 

When a P-wave arrives, we can observe a maximum 
asymmetry of distribution of the vector magnitude of the 
three-component seismograms. As a result of this asym- 
metry, maximum values of kurtosis and skewness are ex- 
pected to occur when a P-wave arrives and this informa- 
tion can be extracted and used for the P-arrival detector. 
Thus, one essential hypothesis of this study is that the 
normalized vector magnitude will show very high kurtosis 
and skewness magnitude when a P-wave arrives. First we 
perform window-by-window normalization of the vector 
magnitude to reduce the huge variations in magnitude re- 
sulting from the ground motion and obtain a zero-mean 
normalized vector magnitude for each window. The nor- 
malized vector magnitude is then used as input to the 
detection and picking system. The automatic detection 
and picking system makes use of the skewness and kur- 
tosis of the input normalized vector magnitude in order 
to detect the P-arrival. For a sliding time window of 
the seismic vector magnitude the values of kurtosis and 
skewness are calculated and the time at which these val- 
ues become maximum/minimum is used to estimate the 
time of arrival of the P-wave. 

Another important contribution of this study is that 
unlike some previous studies which are applied only on 



the vertical component seismogram and used the deriva- 
tives within the kurtosis and skewness values for their 
correction [3|, our technique uses the ratios within the 
kurtosis and skewness values rather than the derivative 
to make a correction in the P-onset time. This new tech- 
nique supposes that the time for P-onset occurs when the 
HOS (skewness and kurtosis) values of the vector mag- 
nitude begins to change drastically during the course of 
gaining their maximum values. In the next subsections 
we discuss the mathematical background of the imple- 
mented technique, our proposed methods, results, dis- 
cussion, and some concluding remarks. 
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FIG. 2: Angle of incidence i for an incoming P-wave approach- 
ing the seimic sensor of a 3-component seismic station. 



II. RATIONALE AND MATHEMATICAL 
BACKGROUND 

A. Background for Vector Magnitude and 
Normalization 

Broadband seismograms represent either the ground 
velocity or displacement. If the three component seismo- 
grams representing the motion of a seismic sensor in the 
three perpendicular directions have magnitudes x, y and 
z, where x and y are usually the east- west and north- 
south components, and z is the vertical component, the 
magnitude of the resultant velocity/displacement vector 
V (Figure 1) can be given by 



V = \f. 



(1) 
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FIG. 1: Three component seismograms (top three panels) 
and their vector magnitude v(n) (the bottom panel) as an 
example. 

For a sequence of digital seismograms, v can be calcu- 
lated sample-by-sample as 



v(n) — \/ x^in) + ipin) + z^(n). 



(2) 



where n is the nth sample. The angle of incidence for 
the seismic ray in Figure 2 can be given in terms of x, y, 
and z component values as 



tan 



\fx^ 



(3) 



where \J x^ -I- y^ represents the total horizontal compo- 
nent of the seismic P-wave velocity. For teleseismic earth- 
quakes (events), the P-wave arrives the seismic sensor 
nearly vertically and thus this angle of incidence is very 
small. For local and regional earthquakes this angle of in- 
cidence becomes significant. The method developed here 
is an improved version of the P Arrival Identification- 
Skewness/Kurtosis technique^- That method is based 
on the observation that the noise in the vertical com- 
ponent seismogram shows zero-mean Gaussian behavior 
until the P-wave arrives. On the other hand, unlike the 
single components, the seismic vector magnitude, v(n) is 
always greater than or equal to zero, and it may not show 
zero mean Gaussian behavior. But we can normalize this 
vector magnitude on a window-by-window basis not only 
to have a zero mean Gaussian behavior, but also to bet- 
ter look at the differences in the computed skewness and 
kurtosis. The normalized vector magnitude is obtained 
through the following transformation: 



Vi(n) 



Vi{n) 



(4) 



Each normalized variable Vi{n) is a rescaled (trans- 
formed) variable of sample Vi{n) for an ith window of 
mean value Vi and standard deviation a.i. This trans- 
formation is linear [5]. Figure 3 displays the Gaussian 
behavior of the normalized vector magnitude for a 3- 
component seismogram. It also shows how the distri- 
bution and skewness and kurtosis values for the normal- 
ized vector magnitude change drastically when P arrives 
(Figure 3(c)). As Figure 3 clearly displays, skewness and 
kurtosis receive much higher values when a P-wave arrival 
is included in their calculation (Figure 3(c)) as compared 
to their near-zero values when P-arrival is not included 
in the shding window (Figure 3(a),(b), and (d)). 
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FIG. 3: A histogram of normalized vector magnitude of 3- 
component seismograms (for 1994 Rukwa seismic event of 
Tanzania recorded by temporary seismic station Hale in Tan- 
zania, East Africa). A window size of 500 samples is used here, 
(a) and (b) are for pure background noise, window starting at 
sample 100 and 150 samples, (c) is when P-wave arrives, and 
so it includes a lot of background noise and P-arrival signal 
around the end of the 500 samples window starting at 931 
sample, and (d) is 5000 samples after P-arrival. The skew- 
ness and kurtosis values for these four windows are as follows: 
(a) skewness — -0.1685, kurtosis — -0.0987 (b) skewness — 
-0.2278, kurtosis = -0.1815 (c) skewness = -8.7562, kurtosis 
= 93.8140 (d) skewness = 0.0350, kurtosis = -0.5883. 



Figure 4 depicts clearly that the kurtosis and skew- 
ness values take their peaks just after P-arrives. This 
normalization process seems to lead to enhancement of 
differences and sets the limit in the range of variations. 
Figure 5 shows the enhancement in the maximum val- 
ues of the quantities under investigation when we apply 
normalization as compared to the values without normal- 
ization. 



B. Mathematical Background for Higher Order 
Statistics 



Skewness (sk) and kurtosis (ku) for a finite length se- 
quence v(n) are estimated by: 
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FIG. 4: (a) Skewness and (b) Kurtosis values for the nor- 
mallized vector magnitude of 3-component seismograms. 
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FIG. 5: Skewness and Kurtosis of normalized and unnormal- 
ized vector magnitude of 3-component seismograms. All hor- 
izontal axes represent sample number. 



Here rhs and as are the mean and standard deviation 
estimates of v{n) and K is the length of the finite se- 
quence. Generally, HOS values for a Gaussian distri- 
bution are zero. In this case, it is demonstrated that 
the normalized magnitude yields higher HOS values for 
P-arrival than the values calculated for the background 
noise and P-wave coda (Figure 3(c)). 



III. HYPOTHESES AND PROCEDURE 

A generalized hypothesis applicable to local and re- 
gional events is proposed. This hypothesis is based on 



the combination of all three-component seismograms in 
a single station, and it is supposed to outperform those 
techniques that are based on just a single component seis- 
mogram. Following the mathematical equations in sec- 
tion 2, the vector magnitude of three component seismo- 
grams for a single station was calculated. This is followed 
by normalization of N-sample window of the vector mag- 
nitude (Eq. 4). While this window slides to the right by 
one-sample at a time the window next to it is automati- 
cally normalized. This normalization continues until the 
last window is normalized. The skewness and kurtosis of 
each normalized window were computed. 

Pure background noise and the seismic signal away 
from P-phase arrival follow a Gaussian distribution with 
nearly zero HOS values. Because of high asymmetry and 
non-Gaussian distribution introduced by the P-arrival 
on the Gaussian background noise, the window with P- 
arrival follows a highly skewed distribution. Maximum 
(Minimum) HOS values are attained for the window that 
includes the arrival of the P-onset and just few additional 
samples of the P-wave coda. Thus, to determine the P- 
onset time more accurately, a correction scheme is intro- 
duced for the additional samples included after the actual 
P-onset. 

The PAI-S/K method of [3j], using the vertical com- 
ponent alone, proposes to use the location of the maxi- 
mum slope as the P-onset time. In this study, however, 
the P-onset time is taken as the time when the HOS 
values of the normalized vector magnitude start to in- 
crease sharply. The HOS values change drastically when 
P breaks out from the background noise level during the 
P-onset. Thus, the magnitude of the ratio of the skew- 
ness value on the right-hand-side at the P-onset to the 
skewness value on the left-hand-side (background noise) 
should be the maximum of all the ratios within the win- 
dow. The same is true with the ratio of kurtosis values 
(Figures 6 to 9). 
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FIG. 6: Application of maximum derivative of skewness for 
correcting P-arrival time according to a 2002 study; it pro- 
vides a correction of 5 samples for the P-arrival. 
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FIG. 7: According to this new study, RHOS value of skewness 
gives a correction of 11 samples for the P-arrival, which is 
by far better than the 2002 study (Figure 6). The correct 
P-arrival happens 12 samples away from the sample of the 
maximum skewness, as shown in Figure 10. 



IV. RESULTS AND DISCUSSION 



A. Experimental Setup 

For this particular study, we used seismic data sets 
obtained from Integrated Research Institutions for Seis- 
mology (IRIS) Data Management Center (DMC) de- 
pository website. Events from several geological ter- 
rains were selected. The seismic networks contributing 
to the IRIS/DMC which are utilized in this study in- 
clude the IRIS-PASSCAL Tanzania Broadband Experi- 
ment (1994/1995), the Global Seismic Network (GSN), 
IRIS/USGS network, PFO of the IRIS/IDA, Southern 
California Seismic Network (SCSN), and Pacific North- 
west Seismograph Network (PNSN). 



B. Results of the new scheme 

Events within regional and local distance range of 50 
km to 1400 km were selected. Only broadband seismic 
data are used in this study. Seismic records with different 
levels of noise were included. The signal analysis was 
performed using Matlab 7.6, R2008a (3 GHz Dual-Core 
Intel mac pro). The performance of the method using 
the seismogram vector magnitude has been compared to 
the PAI-S/K (Figures 6 to 9). The improvement of the 
arrival detection when applying the ratios of HOS values 
as compared to using just the HOS values directly as 
in [3] is clearly shown for both skewness and kurtosis. 
Further discussion is given in the next subsection. Sizes 
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Using maximum derivative in kurtosis, according to 
study, gives a correction of 1 sample for the P-arrival. 
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FIG. 9: According to this new study, applying RHOS 
value of kurtosis offers a correction of 11 samples for the P- 
arrival, which is consistent with the skewness arrival coorec- 
tion scheme of Figure 7, and a much better correction than the 
2002 correction (Figure 8); Moreover, the corrections based 
on maximum derivatives of kurtosis and skewness in the 2002 
study do not agree with each other (Figures 6 and 8), while 
the corrections by RHOS values of this study agree very well 
(Figures 7 and 9). 



technique is applicable to all distance ranges, it is more 
advantageous for regional and local distance range stud- 
ies. Another important difference between the PAI-S/K 
and the new technique is the application of the normal- 
ization scheme that helps to reduce the variation in the 
computed values of skewness and kurtosis from window- 
to- window. 

The third difference, which is a striking one, between 
the PAI-S/K and the technique in this paper is a correc- 
tion procedure, procedure for making correction to the P- 
wave arrival time estimate. To the best of our knowledge, 
this new procedure is introduced in this study for the first 
time to solve such problems. The absolute value of the 
ratio of right-hand-side (RHS) to left- hand- side (LHS) 
HOS values gives a very good estimate of the accurate 
arrival correction. Saragiotis et al. [31] PAI-S/K method 
uses a maximum slope correction procedure. Figures 6 
to 9 show a comparison between these two approaches 
and our results indicate clear improvements when using 
the absolute values of the ratios of the adjacent values of 
HOS (abs(RHS/LHS) of HOS values) as compared to us- 
ing the slope (derivative) of the HOS values. On Figures 
6 and 7, the new ratio method for skewness implemented 
here suggests 11 sample correction while the maximum 
slope in PAI-S/K method suggests 5 sample correction. 
On Figures 8 and 9, the ratio approach for kurtosis sug- 
gests 11 sample correction while the slope approach of 
PAI-S/K technique indicate 1 sample correction. The 
actual correction required can be seen by closely looking 
at the P-arrival very closely (Figure 10). 

Figure 10 indicates that the required correction is 13 
samples. Our approach gives not only a much better cor- 
rection in both kurtosis and skewness cases, but also both 
skewness and kurtosis give consistent correction values. 
Like PAI-S/K we suggest to apply both skewness and 
kurtosis together to detect a P phase arrival. The use 
of both these statistical quantities instead of just one for 
detection will constrain and give a better result than us- 
ing just one of these quantities. The number of false 
alarms also decreases with the use of the two quantities 
simultaneously than using just one of them. 



CONCLUSION 



of sliding window ranging from 30 to 500 samples are 
found to exhibit best results. 



C. Discussion 

Comparison with PAI-S/K and other schemes has been 
made. There are three important differences between the 
new technique presented in this paper and PAI-S/K Q 
method. The PAI-S/K method is among the many ver- 
tical component utilizing methods while our method is 
a three-component based technique. Though the new 



This article has attempted to make use of vector mag- 
nitude of three component seismograms in order to im- 
prove P-wave arrival detection system. Many single com- 
ponent (usually the vertical component), single-station 
based methods have been developed for a P-wave ar- 
rival detection and picking. Since seismic waves from 
regional events approach a seismic station (sensor) in a 
more horizontal incidence than seismic waves from tele- 
seismic events, the energy is generally distributed among 
the horizontal and vertical recording channels. Thus, it is 
not only advantageous to develop a method that makes 
use of all the three components recorded by the single 
station in a combined form but it also gives the tech- 




FIG. 10: Closer look at the vertical (z) component seismo- 
gram including part of the P-arrival; P-arrived about 12 sam- 
ples ahead of the end of the displayed seismogram. 



nique more generality. Though the method is applicable 
to all distance ranges, it is more advantageous for re- 
gional and local distance range studies. We investigated 
the application of the normalization to vector magnitude. 
In contrast, the HOS values for a P-phase arrival do not 
improve as compared to other possible seismic phase ar- 
rivals which may suggest that normalization may play an 
important role in the identification of the-more-difficult- 
to-identify smaller seismic phases. 

We have proposed to apply kurtosis and skewness on 
the combination of three component seismograms for 
picking P wave arrivals automatically and determined 
that the method introduced by Saragiotis et al. [3J can 
be improved if we use the vector magnitude of seismo- 
grams for regional events. This technique makes use of 
the information from all the three component data as 



compared to the application of the technique on the sin- 
gle component, vertical, seismogram. The three compo- 
nent seismograms help us determine the magnitude of 
the total ground motion through their vector magnitude. 
A new approach for making correction of P-wave arrival 
time is also proposed and implemented. The absolute 
value of the ratio of right-hand-side (RHS) to left-hand- 
side (LHS) HOS values gives a very good estimate of the 
accurate arrival correction. This has been compared to 
the derivative (slope) values correction used in the Sara- 
giotis et al. Ql PAI-S/K method and our results indicate 
improvements when using abs(RHS/LHS) of HOS values 
as compared to using the slope of the HOS values. 

We also propose to use both skewness and kurtosis to- 
gether to detect a seismic phase arrival. The use of both 
these statistical quantities instead of just one for detec- 
tion will constrain and give a better result than using just 
one of these quantities. The number of false alarms also 
decreases with the use of thresholds on both quantities 
simultaneously than using just one of these two quanti- 
ties. 

This study has indicated that the RHOS technique can 
be applied on vector magnitude of three component seis- 
mograms for improved detection of P-wave arrivals. It is 
also shown that the new ratio of HOS values technique 
can be applied just on the vertical component seismo- 
grams in lieu of the vector magnitude of three component 
seismograms and it still gives an improved detection of P- 
wave arrivals as compared to that of PAI-S/K. This new 
procedure also provides better detection and picking of 
P-wave arrivals which is essential for locating earthquake 
sources more accurately. The source of the earthquake 
is the point where slippage between the fault surfaces or 
faulting starts inside the earth. Thus, RHOS would en- 
able us to make more accurate earthquake location using 
seismic signals. We strongly believe that this new tech- 
nique has the potential to be applied in a similar fashion 
for accurate detection and location of fractures in ma- 
chines or mechanical systems using acoustic signals. 
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